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Molecular dynamic simulations are performed to investigate a long-time evolution of different 
type initial signals in nonlinear acoustic chains with realistic Exp-6 potential and with power ones. 
Finite number of long-lifetime kink-shaped excitations is observed in the system in thermodynamic 
equilibrium. Dynamical equilibrium between the processes of their growth and decay is found. 
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I. INTRODUCTION 

Power expansions of interatomic potentials (such as 2-4 
and 2-3-4 potentials) are frequently used for investigation 
of strongly nonlinear phenomena ||, |[ ^, ||, ^ . How- 
ever, for soliton-like excitations with large amplitudes 
one can not be sure that high-order anharmonisms may 
be omitted. So, realistic (e.g. Lennard- Jones, Morse) 
potentials were used to take all-order anharmonisms into 
account [||, |. 

Several types of soliton-like excitations were found in 
discrete monoatomic chains with anharmonic intersite 
potentials (acoustic chains). 

Breathers - time-periodic space-localized one- 
parameterl^l modes with frequency as a parameter 
are revealed in FPU chains with 2-4 potential p^ . 
Another solutions - dark solitons (phase shift in stable 
time-periodic extended nonlinear mode) are proved to 
exist in FPU chains jl^, In monoatomic acoustic 

chains with realistic Lennad- Jones and Morse potentials 
breathers are not observed in numerical experiments, 
and, moreover, they are forbidden However, they 
arise if time-periodic driving field is applied to the chain 



Moving kinks - time non-periodic solitons with step- 
shaped profiles of atom displacements |l^, 15, 16, 17] 
and pulses - solitons with bell-shaped profiles |15, 18, 1£] 
are predicted in acoustic chains. These solitons are one- 
parameter ones with a velocity as a parameter. Numer- 
ical results for pulses were obtained in g |l^ for FPU 
and Toda chains. Exact non-topological moving kinks 
are well known in integrable acoustic Toda chain [p^ . 
Existence theorems of non-topological mo ving k inks are 

and for 



also proved for non-integrable FPU chains 1 14 



acoustic chains with interatomic potentials of arbitrary 
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[1] Necessary conditions for existence of two-parameter movini 
breathers (frequency and velocity are a parameters) and numer- 
ical algorithm for their searching are proposed in H] 



powers We don't know other solutions for moving 
kinks in discrete case for acoustic chains. In continual 
(long wave) limit moving kinks are observed in real ul- 
trasonic experiment and investigated numerically taking 
third-order anharmonism into account (KDV equation) 

Thus, no any soliton-like solutions have been obtained 
in discrete case for non-integrable acoustic chains with re- 
alistic interatomic potentials. In this paper spontaneous 
creation of discrete moving kinks is observed in numerical 
experiment for the chain with such potential. 

The aim of this article is molecular dynamics study of 
discrete acoustic chain with realistic interatomic poten- 
tial at high temperatures after thermalization. The role 
of different orders of anharmonicity in power expansion 
of realistic potential is investigated. 

The methods used are described in the Sec. II. Numeri- 
cal results and discussions for thermalized gas of sponta- 
neously appearing moving kinks are given in the Sec. HI. 



II. METHODS AND APPROXIMATIONS 



Equations of motion 
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are integrated numerically with periodic boundary con- 
ditions for acoustic chain containing 100 atoms. Nearest 
neighbor approximation is used. Analogous results are 
obtained when six neighbors and 1500 atoms are taken 
into account. 

We use the system of units in which energy is mea- 
sured in [K], mass in [a.m.u.] does, and distance is mea- 
sured in the units of equilibrium interatomic distance 
[do]. The unit of time is equal to the period of har- 
monic zone boundary phonon mode (Tq — 8.3716 • 10"^'^ 
sec) . It corresponds to the case of neon with atom mass 
771=20.18 a.m.u. Velocities are measured in the units of 
sound velocity Vs = 12.72[{K /a.m.u.y^^]. 
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Six-order symplectic Yoshida algorithm [gO| is applied 
for numerical simulations. Time step is equal to 0.008 
To which provides energy conservation with the accuracy 
l\E/E = 10"^ during whole simulation time {t = IO^Tq). 

Realistic interatomic Exp-6 potential 

y-p-6(|^^ - r,\) = A, exp(-a(a;., - 1)) - ^.t'-^; 

Xij = |ri - Tjl/do (2) 

is used with the parameters Ao=35.9335 K, a=13.6519 
obtained ab — initio for the dimer of neon [ pT| . Single 
empirical parameter do — 3.091A is fitted to equilibrium 
interatomic distance in neon dimer p^ . 
Four types of initial conditions are used. 

• Zone boundary mode with wave vector fc = vr (tt- 
mode). Equal opposite initial velocities |Vo| = 
Q.2Avs were assigned to neighbor atoms. 
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and taking higher powers into account consequently. 



III. RESULTS AND DISCUSSIONS 

Excluding the cases specially mentioned, acoustic 
chain with Exp-6 potential is considered with different 
initial conditions in all the subsections. Chains with 
power and Toda potentials are also reported. 



A. TT-mode 

To the time w SOTq initial 7r-mode perturbed by nu- 
merical noise is destroyed by period-doubling instability 
(see also Q). In the following, the motion of atoms be- 
come more and more chaotic. Simultaneously, the pro- 
cess of space concentration of energy takes place. Well- 
localized solitons with the width compatible to lattice 
constant are spontaneously created since ~ 6OOT0. Thus, 
the solitons in the chain with realistic potential appear 
after thermalization. On the contrary, breathers in FPU 
chain arise directly as a result of 7r-mode destruction, and 
the system is termalized after decay of breathers. 

The thermalization is defined here as a state when the 
energy distribution of atoms obtained in numerical simu- 
lation is close to Gibbs one. The system persists in such 
state for a long time. The deviations from Gibbs distri- 
bution caused by solitons are possible at high-energy tile 
only. 



• White noise. Random initial velocities were fixed. 

• Shock waves. Initial velocities |Vb| = 0.96ws di- 
rected inside the chain were given to three atoms 
at each end of the chain. 

• Exact breathers. Initial displacements of atoms 
were fixed according to breathers exact form 
A(. . . , 0, -1/2, 1,1/2, 0, . . .) for odd-parity breather 
and A(. . . , 0, — 1, 1, 0, . . .) for even-parity one [||. 
The amphtude is A=0.279 do- 

In all the cases the same energy _E=90 K/atom is fed 
into the system by the initial conditions. This energy is 
two times larger then cohesive energy of neon dimer. 

To investigate the role of different order anhar- 
monisms, the simulations have performed on acoustic 
chains with power interatomic potentials. These poten- 
tials are defined expanding Exp-6 potential 



(3) 
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In the Fig. g energy distribution of atoms E{£) = en{e) 
is given. The E{e) is defined so that dE{e) = E{e)de 
is total energy of atoms which fall into energy interval 
(e, e + de). The n{e)de is the number of atoms with the 
energies residing in the same interval. Solid curve E{e) 
corresponds to Gibbs distribution n(e) obtained in har- 
monic approximation. The circles denote the distribution 
E{e) observed in numerical experiment, when one of the 
most intensive solitons appears. This soliton gives a local 
peak on the tile of experimental E{e). The position of 
the peak on energy scale manifests that the energy con- 
centrated in the soliton is « 15 times as mush as average 
energy of an atom in the system. At all the rest ener- 
gies experimental distribution agrees with Gibbs one. It 
indicates that majority of atoms move chaotically. 

Fractional deviation of experimental distribution n{e) 
from Gibbs one is reported in Fig. ^ (numbered by 1). 
Experimental and Gibbs distributions E{e) are also given 
here (numbered by 2). The value of experimental n(e) is 
w 36 times as much as Gibbs one near the peak at the 
tile, while Unexper / noib ~ 1 at all the rest energies. 

Typical pattern of map of tracks of solitons is plotted in 
Fig. y where single powerful soliton is seen clearly. Time 
periodicity in the motion of atoms is not observed. Ve- 
locity of the soliton (the slope of the curve of energy max- 
imum position via the time) is equal to v — 2.6vs. Nu- 
merical experiment shows that the the velocities and the 
amplitudes of observed solitons are connected uniquely. 
Therefore, these solitons are one-parameter ones. Be- 
sides, they are pass each through others without loss of 



y-p-6(|^^ - r,|) = Vo + -f{x,, If + -f{x., - 1)^ + -f{x,, - if + . . . 



3 




S, K 

FIG. 1: Energy distribution of atoms E{e) — e7i{e). Solid curve corresponds to Gibbs form of the n{e), the circles denote the 
distribution E{e) obtained in numerical experiment. To get experimental distribution the energies of atoms are averaged for 
the time 11 To. The distribution corresponds to one of the most intensive kinks {t = 188947ro). 



their individuality like to solitons in integrable systems 
(Fig. I). 

The form of obtained solitons can be determined from 
enlarged fragment of the map of tracks plotted in Fig. |^. 
The velocities of atoms are given here. The atom in the 
soliton moves in the same direction that the soliton does 
and transfer its momentum to the next neighbor. It can 
be said with good accuracy that this soliton is localized 
on two sites only. Space dependence of atoms velocity 
has bell-shaped form at each moment of time. Therefore, 
the displacements of atoms have step-shaped form, and 
the solitons can be identified as non-topological moving 
kinks. 

Let use the parameter 

Co^NY,Ef/{Y.E,r (4) 

i=l 1=1 

as quantitative indicator of energy concentration H| . The 



Cq—I if all the atoms have the same energies, Co=1.75 
for equilibrium Gibbs distribution in harmonic approxi- 
mation [|, and Cq= N if all the energy is concentrated 
on one site only. Parameter Co indicates energy con- 
centration in all the solitons existing simultaneously, but 
speaks nothing about the number of solitons and their 
individual contributions. In Fig. ^ the Cq via the time 
is given for the chain after thermalization. For the time 
interval under consideration Co oscillates around its av- 
erage value in the frames 1.7-3.3, and it doesnt exceed 
1.5-4.2 during whole simulation time {t = IO^Tq). It 
follows from this oscillating behavior of the Cq that the 
probabilities of kinks growing and decay are equal. 

Different periodic components of Co{t) are separated 
by fast Fourier transform which is performed taking 
8192 time points into account. To suppress weak high- 
frequency fluctuations of the Co caused by thermal mo- 
tion of atoms initial time dependence of Co was smoothed 
by adjacent averaging using 25 time points. Amplitude 
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FIG. 2: Fractional deviation of experimental distribution n{e) 
from Gibbs one, nexper{£)/noib{£), ai t — 2316248To (1 - 
black sircles connected by the line; right scale). The digit 
2 denotes experimental and Gibbs distributions E{e) at the 
same t. 



Fourier spectrum smoothed using 5 frequency points is 
plotted in Fig. ^ (the curve 1). A number of considerable 
peaks is observed at the frequencies 10~^ — 10~'^(To)~^. 
However, the width of these peaks is compatible to the 
distance between them, and each peak by itself is seen 
to be meaningless. Thus, one can evaluate maximal fre- 
quency (minimal period) of long-time oscillations of Cq (t) 
only. For convenience we smooth frequency dependence 
of Fourier amplitudes using more number of points (50 
points, the curve 2 in the Fig. We define maximal fre- 
quency as a frequency at which Fourier amplitudes drop 
two times. Then, minimal period of Co{t) oscillations 
is evaluated to be ~ IOOOTq, which agrees with the life- 
time of most intensive kinks defined directly using tracks 
maps. 

Numerical simulations show a number of kinks to ex- 
ist simultaneously in the chain. Some of them persist on 
the stage of growing and others are on the stage of decay. 
Moreover, there are a long-time fluctuations of number of 




FIG. 3: Map of tracks of moving kink from numerical exper- 
iment of the Fig. 0. Horizontal axis indicates the position 
along the chain, vertical axis corresponds to the time (time 
is going downward). The energies and velocities of atoms are 
denoted as short vertical lines with the length proportional to 
the magnitude (positive direction is downward). 



kinks which correlate with the fluctuations of energy con- 
centration parameter Cq: the more the magnitude of Cq 
the less the number of kinks^^l . For example, for strong 
concentration of energy (Co=3.7) corresponding to the 
Fig. ^ single kink runs along the chain. The number 
of kinks remains, practically, unchanged when enlarged 
amounts of energy of 300 K/atom and 1000 K/atom are 
initially fed into the chain, while the energy concentrated 
in each kink enhances. 

In conclusion, there is an hierarchy of excitations in 
the chain after thermalization, starting with most strong 
well-localized and long-lifetime kinks and finishing with 



[2] Similar correlation was observed on acoustic chain with Lennard- 
Jones potential . 
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FIG. 4: Interaction of moving kinks. 
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FIG. 5: The form of moving kink from the Fig. |3|. Time step 
is equal 0.022 Tq. 



phonons. Weak excitations deserve as the background 
with which strong kinks interact. Time after time strong 
kinks appear spontaneously from the background. For 
certain time each kink exhibits systematic tendency to 
growth, which changes by stable tendency to decay. The 
behavior of ensemble of such kinks leads to oscillating 
time dependence of energy concentration parameter, and 
all the properties averaged for a long enough time are 
steady. 



B. White noise 

Although there is a stage of chaotic motion of atoms 
after decay of initial mode the question remains about ef- 
fect of initial conditions on the form of created solitons. 
To exclude this effect we use white noise as initial condi- 
tions. The system will choose most appropriate solitons. 

At the beginning of the simulations atoms move chaot- 
ically, but the distribution E{e) is differ a bit from Gibbs 
one. The process of energy concentration starts imme- 
diately. After fti 80To the E{e) agrees with Gibbs one. 
Well localized kinks are created since « 400To. Their 
form and properties are identical to ones for the kinks 
arising after 7r-mode decay. 



Thus, whatever initial conditions, acoustic chain with 
realistic Exp-6 potential comes to the same state char- 
acterized by finite number of long-lifetime high-energy 
supersonic moving kinks. 

C. Shock waves 

The form of spontaneously appearing moving kinks is 
similar to the form of shock waves. We adopt shock waves 
as initial conditions, which permit to set the amplitudes 
of the kinks in controlled manner. 

Instantly each initial shock wave including three atoms 
splits into three kink-like excitations moving with the ve- 
locities 4.5-5 Vs', Co=17.4. They are locaHzed on two sites 
only. Kink-like excitations lose the energy slowly. To 
w llOOTo they become slightly delocalized in space: two, 
three or more atoms moving in the same direction are 
contained in the kinks. The parameter Cq drops to 3.1, 
kinks velocities decrease to 2.3-2.8 Vs- So, quantitative 
characteristics of shock waves become compatible to ones 
for spontaneously created kinks. 

Numerical experiments on shock waves with different 
initial amplitudes show that the more energy is stored in 
the wave the more its lifetime. 

We consider the propagation of the same shock waves 
in integrable Toda chain where non-topological moving 
kinks are exact solutions. Like FPU chain, each initial 
shock wave splits into three kinks with the form, veloc- 
ities and parameter Co similar to ones for Exp-6 poten- 
tial. In the following these kinks transfer a bit of energy 
to background atoms, and, further, neither the form nor 
the energies of the kinks are altered to the end of simu- 
lation (t=28000 To). 

D. Exact breathers 

It was found that in monoatomic acoustic chains with 
realistic Lennad-Jones and Morse potentials breathers 
are not created We have performed numerical sim- 
ulations on the chain with Exp-6 potential with initial 
displacements of atoms in the form of single odd-parity 
breather placed in the center of the chain. For the time 
less then To the breather decays into two kinks with equal 
energies moving in opposite directions. In the following 
shock wave scenario is realized. Kink-like excitations are 
delocalized slightly in space, and to 5000 To the parame- 
ter Co and kinks velocities decrease and become compat- 
ible to the characteristics of spontaneously created kinks. 
Further, each kink grows and decays as it was described 
above. The simulations with even parity breather give 
analogous result. 

Instability of breathers in the chain with Exp-6 poten- 
tial is caused by strong asymmetry of the potential. In- 
side the breather, one of interatomic bonds is compressed 
and its potential energy is increased. If the potential 
would be symmetrical, increasing of the potential energy 
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FIG. 6: Parameter of energy concentration via the time (black line). White line is average value < Co >. 
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E. Power potentials 

The form and some properties of moving kinks in the 
chain with Exp-6 potential are similar to ones of the so- 
lutions of Bussinesk equation. It can be obtained from 
Eq. in continual approximation if cubic anharmonism 
only is taken into account. This similarity would suggest 
that high order anharmonisms could be neglected. To 
clarify the role of different order anharmonisms we per- 
form the simulations for the chains with the potentials 
obtained as power expansion Eq. ^ of Exp-6 potential. 

The next cases are considered (initial conditions - tt- 
mode). 
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FIG. 7: Amplitude Fourier spectrum of Co(t). The curve 1 
is the spectrum smoothed using 5 points; the curve 2 is the 
spectrum smoothed using 50 points. 



would be equal when the atoms both approach each other 
and remove. So, if the slope of potential curve is large 
enough, the pair of atoms would vibrates almost like iso- 
lated bound state as it does in FPU chain. In Exp-6 
potential attractive part is much weaker then repulsive 
one, and the atoms, if approached, remove almost freely. 
They collide with neighbor atoms and transfer whole mo- 
mentum to them. As a result, two kinks run along the 
chain in opposite directions. 



0,005 • 2-3 potential. 7r-mode decays, the solitons appear 
and collapse immediately. 

• 2-4 and 2-3-4 potentials. 7r-mode decays by mod- 
ulational instability, and the breathers occur in 
agreement with [0, g, ^, ^, ||, ^ . Maximal velocity 
of breather is ~ O.lvs, the lifetime is w 120To (2-4 
potential) and w 50To (2-3-4 potential). Unlike the 
kinks, breathers interact each with other strongly. 
After breathers decay movement of atoms becomes 
chaotic, energy distribution agrees with Gibbs one. 
In the following kink-like solitons are spontaneously 
created in the chain^'^l . However, lifetime of these 
kinks is very short (« 50ro)i indicating that both 
kinks and breathers interact strongly in the chains 
with 2-4 and 2-3-4 potentials. 



[3] Numerical simulation shows kink-like solitons to arise in the 
chains with 2-4 and 2-3-4 potentials regardless of initial con- 
ditions. 
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• 2-3-4-5-6 and 2-3-4-5-6-7-8 potentials. 7r-mode is 
stable. 

• 2-3-4-5-6-7-8-9-10 potential. Tr-mode decays trough 
period-doubling instability. Solitons identical to 
kinks in the chain with Exp-6 potential are created 
(with long lifetime and weakly interacted each with 
others). 

Thus, in the chain with 2-3 potential nonlinear term is 
not balanced by dispersion one (at least at the energy 90 
K/atom), and solitons are collapse. Cubic anharmonism 
added to 2-4 potential give only quantitative changes in 
the picture of solitons formation. The 5-6 and 7-8 powers 
stabilize 7r-mode. Higher (9-10) powers give the change 
of the type of instability, and the properties of kinks coin- 
cide with ones in realistic Exp-6 potential. So, quantita- 
tively small corrections give qualitative change of atoms 
behavior, and realistic interatomic potentials should be 
used to avoid such difficulties in real crystals. 

IV. SUMMARY 

The behavior of nonlinear acoustic chain with realis- 
tic Exp-6 potential is investigated. Regardless of the type 
of initial conditions, moving non-topological kink-shaped 
solitons occur spontaneously after thermalization of the 
system. These kinks are one-parameter solitons with the 
velocity as a parameter. They interact each with oth- 
ers weakly and pass each through others without loss 
of their individuality like solitons in integrable systems. 
Traveling along the chain, each kink, firstly, collect the 



energy from the background atoms, and, in the following, 
transfer it to small-amplitude phonons. Dynamical equi- 
librium between the processes of kinks growth and decay 
takes place. So, there is finite number of high-energy ex- 
citations in the system at each moment of time. Time 
intervals when the most amount of energy is concentrated 
in the kinks are repeated periodically. Average lifetime 
of a kink is estimated to be ~ IOGOTq. 

Acoustic chains with power interatomic potentials are 
investigated adding higher powers of atom displacements 
consequently up to ten powers. In the chains with 2-4 
and 2-3-4 potentials discrete breathers are observed in 
agreement with j^, |^, ^, ^ |5[ |^. After breathers decay 
thermal equilibrium is achieved and kink-like solitons also 
appear in these chains. However, these kinks interact 
strongly each with others and have short lifetime (50-120 
To). Thus, unlike breathers resided to power potentials 
only, the kinks can be created both in power and realistic 
potentials whatever initial conditions. 

Molecular dynamics study of 7r-mode stability in the 
chains with power potentials shows the change of the type 
of instability (from modulational to period-doubling) 
when high (9 and 10) powers are added. It indicates 
that extreme caution is required using power expansions 
of interatomic potentials in real crystals. 

Thus, in thermal equilibrium long-lifetime high-energy 
supersonic moving kinks appear in acoustic chain with 
Exp-6 potential at high temperatures. One can expect 
that similar excitations can exist in three-dimensional 
case and give considerable enhancing of thermal conduc- 
tivity at high temperatures and also effect on long-time 
strength of solids. 
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